library(ggrepel)
library(clusterProfiler)
library(ggplot2)
library(enrichplot)
library(readxl)
library(plotly)
library(VennDiagram)

library(RColorBrewer)

source('pathway_function.R')
color = colorRampPalette(rev(brewer.pal(n = 7, name =
                                          "RdYlBu")))(100)

genelist=c('TP53', 'BAX', 'BAK1', 'OTUD5', 'PMAIP1', 'FLT3', 'CDK2', 'UBE2J2', 'MFN2', 'TOMM70A', 'MARCH5')

files=c('./data/crispr/Molm13_MCL1i_DMSO_Day16.gene_summary.txt',
        './data/crispr/Molm13_combo_DMSO_Day16.gene_summary.txt',
        './data/crispr/VA_DMSO_Day16.gene_summary.txt',
        './data/crispr/Ven_DMSO_Day16.gene_summary.txt',
        './data/crispr/Aza_DMSO_Day16.gene_summary.txt')
axis_names=matrix(c('MCL1i lfc','MCL1i & Ven lfc','MCL1i lfc','Ven & Aza lfc','MCL1i & Ven lfc','Ven & Aza lfc',
                    'Ven lfc','MCL1i & Ven lfc'),2,4)
comparisons=matrix(c(1,2,1,3,2,3,4,2),2,4)
totdata=list()

for (kk in 1:length(files)) {
  temp=read.csv(files[kk],header=T,sep='\t',stringsAsFactors = F,check.names = T)
  totdata[[kk]]=temp
}
for (k in 1:ncol(comparisons)) {
  d1=totdata[[comparisons[1,k]]]
  d2=totdata[[comparisons[2,k]]]


  d1$id[d1$id=='5-Mar']='MARCH5'
  d2$id[d2$id=='5-Mar']='MARCH5'
  genes=intersect(d1$id,d2$id)
  d1=d1[match(genes,d1$id),]
  d2=d2[match(genes,d2$id),]
  for (ii in 1:nrow(d1)) {
    if (d1$neg.lfc[ii]>=0)
      d1[ii,3:8]=d1[ii,9:14]
    
    if (d2$neg.lfc[ii]>=0)
      d2[ii,3:8]=d2[ii,9:14]
  }
  
  genes=union(d1$id[d1[,7]>=0],d2$id[d2[,7]>=0])
  
  d1=d1[match(genes,d1$id),1:8]
  
  d2=d2[match(genes,d2$id),1:8]
  
  plotdata=d1
  
  plotdata[,2]=d1$neg.lfc
  plotdata[,3]=d2$neg.lfc
  plotdata[,2]=scale(plotdata[,2])
  plotdata[,3]=scale(plotdata[,3])
  
  plotdata[,4]=d1$neg.p.value
  plotdata[,5]=d2$neg.p.value
  plotdata[,6]=-log10(sqrt(d1[,4]*d2[,4]))
  plotdata[,7]=0
  plotdata[plotdata[,4]<0.05,7]=1
  plotdata[plotdata[,5]<0.05,7]=2
  plotdata[plotdata[,4]<0.05 & plotdata[,5]<0.05,7]=3
  
  names(plotdata)=c('id','MCL1i_lfc','combo_lfc','Pval1','Pval2','Pval_comb','significance')
  plotdata$significance=factor(plotdata$significance,levels=c(0,1,2,3),labels=c('not sig','x_sig','y_sig','both_sig'))
  plotdata=plotdata[,1:7]
  temp=plotdata
  getlabel <- function(x) {
    for (i in 1:length(ind))
      if (!is.na(as.numeric(x[ind[i]])))
        x[ind[i]]=round(as.numeric(x[ind[i]]),2)
    return(paste('</br>',paste(names(temp)[ind],x[ind],sep=':',collapse='</br>'),sep='',collapse='') )}
  ind=1:6
  plotdata$label=apply(plotdata,1,getlabel)
  
  plotdata[order(plotdata$significance),]
  lims=ceiling(max(abs(plotdata[,c('MCL1i_lfc','combo_lfc')])))
  
  
  alpha=0.8
  psize=0.8
  alpha2=1.0
  psize2=1.8
  p1=ggplot(plotdata, aes(x=MCL1i_lfc, y=combo_lfc )) + 
    geom_point(data=plotdata[plotdata$significance %in% 'not sig',],aes(x=MCL1i_lfc, y=combo_lfc,color=significance),alpha = alpha,size=psize,shape=16)+
    geom_point(data=plotdata[plotdata$significance %in% c('x_sig','y_sig'),],aes(x=MCL1i_lfc, y=combo_lfc,color=significance),alpha = alpha,size=psize,shape=16)+
    geom_point(data=plotdata[plotdata$significance %in% 'both_sig',],aes(x=MCL1i_lfc, y=combo_lfc,color=significance),alpha = 0.5,size=psize,shape=16)+
    geom_point(data=plotdata[plotdata$significance %in% 'not sig' & plotdata$id %in% genelist,],aes(x=MCL1i_lfc, y=combo_lfc,color=significance),alpha = alpha2,size=psize2,shape=16)+
    geom_point(data=plotdata[plotdata$significance %in% c('x_sig','y_sig') & plotdata$id %in% genelist,],aes(x=MCL1i_lfc, y=combo_lfc,color=significance),alpha = alpha2,size=psize2,shape=16)+
    geom_point(data=plotdata[plotdata$significance %in% 'both_sig' & plotdata$id %in% genelist,],aes(x=MCL1i_lfc, y=combo_lfc,color=significance),alpha = 1.0,size=psize2,shape=16)+
    geom_point(data=plotdata[plotdata$significance %in% 'both_sig' & plotdata$id %in% genelist,],aes(x=MCL1i_lfc, y=combo_lfc,color=significance),alpha = 1.0,size=psize2,shape=1,color='black')+
    labs(x=axis_names[1,k],y=axis_names[2,k])+
    theme_bw()+
    theme(axis.line.x = element_line(color="black", size = 1),
          axis.line.y = element_line(color="black", size = 1))+
    scale_x_continuous(limits = c(-lims,lims),breaks = seq(-lims, lims, by = 2),minor_breaks = seq(-lims, lims, by = 2))+
    scale_y_continuous(limits = c(-lims,lims),breaks = seq(-2*lims,2*lims, by = 2),minor_breaks = seq(-2*lims,2*lims, by = 2))+
    geom_hline(yintercept=0, linetype="dashed", color = "grey")+
    geom_vline(xintercept=0, linetype="dashed", color = "grey")+
    scale_color_manual(breaks=c('both_sig','y_sig','x_sig','not sig'),values=c('#d63638','#2a9d8f','#edc46a','grey'))+    
    geom_text_repel(data = plotdata[plotdata$id %in% genelist,], aes(label = id),segment.alpha=0,
                    color = "black", size = 4, point.padding = 0.2, min.segment.length = unit(0, 'lines'))+
    theme(axis.text = element_text(size = 15),
          axis.title = element_text(size = 15),
          text = element_text(size = 10),
          legend.text=element_text(size=15),
          legend.title=element_text(size=20),
          panel.border = element_rect(colour = "black", fill=NA, size=2))
  
  if (!dir.exists('./figures'))
    dir.create('./figures',recursive=T)
  
  ggsave(paste('./figures/Figure_1C_',gsub(' lfc|&','',axis_names[1,k]),'_vs_',
               gsub(' lfc','',axis_names[2,k]),'_scatter_filter_50_zscore.pdf',sep='',collapse=''),p1,width=9,height=7) 
  kegg_analysis(plotdata$id[plotdata$significance %in% 'both_sig'],paste0(axis_names[1,k], '_vs_',axis_names[2,k]))
  
  ###doing Venndiagram
  venndata=list()
  venndata[[gsub(' lfc','',axis_names[1,k])]]= plotdata[plotdata$significance %in% c('x_sig','both_sig'),1]
  venndata[[gsub(' lfc','',axis_names[2,k])]]= plotdata[plotdata$significance %in% c('y_sig','both_sig'),1]
  
  venn.diagram(
    venndata,
    #category.names = c("Set 1" , "Set 2 " , "Set 3"),
    filename = paste('./figures/sFigure_1BE_',gsub(' lfc|&','',axis_names[1,k]),'_vs_',
                     gsub(' lfc','',axis_names[2,k]),'_venn_diagram.tiff',sep='',collapse=''),
    fill=c('blue','red'),
    output=TRUE
  )
}
###GO pathway analysis on itersection of 3 condition

temp=read.csv('./data/crispr/intersection_genes_MCL1i_MCL1i+ven_Ven+Aza.csv',header=T,stringsAsFactors = F)
genes=temp[,1]
go_analysis(genes,'MCL1i_MCL1i+ven_Ven+Aza')


